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We construct general anisotropic cosmological scenarios governed by an f{R) gravitational sector. 
Focusing then on Kantowski-Sachs geometries in the case of _R"-gravity, and modelling the matter 
content as a perfect fluid, we perform a detailed phase-space analysis. We find that at late times 
the universe can result to a state of accelerating expansion, and additionally, for a particular n- 
range (2 < n < 3) it exhibits phantom behavior. Furthermore, isotropization has been achieved 
independently of the initial anisotropy degree, showing in a natural way why the observable universe 
is so homogeneous and isotropic, without relying on a cosmic no-hair theorem. Moreover, contracting 
solutions have also a large probability to be the late-time states of the universe. Finally, we can also 
obtain the realization of the cosmological bounce and turnaround, as well as of cyclic cosmology. 
These features indicate that anisotropic geometries in modified gravitational frameworks present 
radically different cosmological behaviors comparing to the simple isotropic scenarios. 
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I. INTRODUCTION 

The observable universe is homogeneous and isotropic 
at a great accuracy [l|, leading the large majority of cos- 
mological works to focus on homogeneous and isotropic 
geometries. The explanation of these features, along with 
the horizon problem, was the main reason of the infla- 
tion paradigm construction 0]. Although the last sub- 
ject is well explained, the homogeneity and the isotropy 
problems are, strictly speaking, not fully solved, since 
usually one starts straightaway with a homogeneous 
and isotropic Friedmann-Robertson- Walker (FRW) met- 
ric, and examines the evolution of fluctuations. On the 
other hand, the robust approach on the subject should 
be to start with an arbitrary metric and show that in- 
flation is indeed realized and that the universe evolves 
towards an FRW geometry, in agreement with observa- 
tions. However, the complex structure of such an ap- 
proach allows only for a numerical elaboration [3|, and 
therefore in order to extract analytical solutions many 
authors start with one assumption more, that is inves- 
tigating anisotropic but homogeneous cosmologies. This 
class of geometries was known since a long time ago (4| , 
and it can exhibit very interesting cosmological behavior, 
either in inflationary or in post-inflationary cosmology 
0. Finally, note that such geometries may be relevant 
for the description of the black- hole interior Q . 

The most well-studied homogeneous but anisotropic 
geometries are the Bianchi type (see [3, Q and references 
therein) and the Kantowski-Sachs metrics [MtIIIj either 
in conventional or in higher-dimensional framework. Al- 
though some Bianchi models (for instance the Bianchi 
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IX one) arc more realistic, their complicated phase-space 
behavior led many authors to investigate the simpler but 
very interesting Bianchi I, Bianchi III, and Kantowski- 
Sachs geometries. In these types of metrics one can ana- 
lytically examine the rich behavior and incorporate also 
the matter content |lll - [l5| , obtaining a very good picture 
of homogeneous but anisotropic cosmology. 

On the other hand the observable universe is now 
known to be accelerating [l^ , and this feature led physi- 
cists to follow two directions in order to explain it. The 
first is to introduce the concept of dark energy (see [13] 
and references therein) in the right-hand-side of the field 
equations, which could either be the simple cosmological 
constant or various new exotic ingredients |18l42(]| . The 
second direction is to modify the left-hand-side of the 
field equations, that is to modify the gravitational the- 
ory itself, with the extended gravitational theories known 
as /(i?)-gravity (see 21, 22] and references therein) being 
the most examined case. Such an approach can still be 
in the spirit of General Relativity since the only request 
is that the Hilbert-Einstein action should be generalized 
(replacing the Ricci scalar R by functions of it) asking for 
a gravitational interaction acting, in principle, in differ- 
ent ways at different scales. Such extended theories can 
present very interesting behavior and the corresponding 
cosmologies have been investigated in detail ^23-27] . 

In the present work we are interested in investigating 
homogeneous but anisotropic cosmologies, focusing on 
Kantowski-Sachs type, in a universe governed by f{R)- 
gravity. In particular, we perform a phase-space and sta- 
bility analysis of such a scenario, examining in a sys- 
tematic way the possible cosmological behaviors, focus- 
ing on the late-time stable solutions. Such an approach 
allows us to bypass the high non-linearities of the cos- 
mological equations, which prevent any complete analyt- 
ical treatment, obtaining a (qualitative) description of 
the global dynamics of these models. Furthermore, in 



these asymptotic solutions we calculate various observ- 
able quantities, such are the deceleration parameter, the 
effective (total) equation-of-state parameter, and the var- 
ious density parameters. We stress that the results of 
anisotropic f{R) cosmology are expected to be different 
than the corresponding ones of /(i?)-gravity in isotropic 
geometries, similarly to the differences between isotropic 
[285 ^nd anisotropic [20| considerations in General Rela- 
tivity. Additionally, the results are expected to be dif- 
ferent from anisotropic General Relativity, too. As we 
see, anisotropic f{R) cosmology can be consistent with 
observations. 

The paper is organized as follows: In section |ll] we 
construct the cosmological scenario of anisotropic /(i?)- 
gravity, presenting the kinematical and dynamical vari- 
ables, particularizing on the f{R) — -R" ansatz in the 
case of Kantowski-Sachs geometry. Having extracted the 
cosmological equations, in section IIIII we perform a sys- 
tematic phase-space and stability analysis of the system. 
Thus, in section ITVl we analyze the physical implications 
of the obtained results, and we discuss the cosmological 
behaviors of the scenario at hand. Finally, our results 
are summarized in section FVl 



II. ANISOTROPIC f{R) COSMOLOGY 



In the following, we will focus on the Kantowski-Sachs 
geometry [I^l, since it is the most popular anisotropic 
model, and since all solutions are known analytically in 
the case of General Relativity, even if some particular 
types of matter are coupled to gravity [III Q^. These 
are spatially homogeneous spherically symmetric mod- 
els [30, |31|, with 4 Killing vectors: du; 



cos(/3 d^ — 
sin if cot 'd d^ , sin ip d^ + cos (p cot d d^,dx [S^l . 

Let us now consider relativistic fluid dynamics (see for 
e.g. [33]) in such a geometry. For any given fluid 4- 
velocity vector field u^\ the projection tensor^ 



i^^u — 9 ^J'l^ 



U^j_Uy 



projects into the instantaneous rest-space of a comoving 
observer. It is standard to decompose the first covariant 
into its irreducible parts 



derivative V^w 



^liUu = 



'U^^iiv + Crf_iu + o^^A"' ^ ^A"'' 



where (T^^ is the symmetric and trace-free shear tensor 
(ct^^ = (^[fiv), (^tii^u" = 0, o-^^ = 0), ujf^^ is the antisym- 



metric vorticity tensor (w^i^ 



^[M^l 



^/ii^ 



0) and Up 



is the acceleration vector defined as u^ — u'^'V^u^ (a dot 
denotes the derivative with respect to t). In the above 
expression we have introduced also the volume expansion 
scalar. In particular. 



In this section we present the basic features of 
anisotropic Bianchi I, Bianchi III, and Kantowski-Sachs 
geometries. After presenting the kinematical and dynam- 
ical variables in the first two subsections, we focus on the 
Kantowski-Sachs geometrical background and on the i?"- 
gravity in the last two subsections. 



e 



V^u^ 



(4) 



which defines a length scale £ along the flow lines, describ- 
ing the volume expansion (contraction) behavior of the 
congruence completely, via the standard relation, namely 



e 



se 



(5) 



A. The geometry and kinematical variables 

In order to investigate anisotropic cosmologies, it is 
usual to assume an anisotropic metric of the form [l3, 
Q: 

ds' = -N{tfdt'' + [ei\t)]-''dr'' + [e2''{t)]-^[d^'' + S{^fd^^], 

(1) 
where l/ei^(i) and l/e2^(t) are the expansion scale fac- 
tors. The frame vectors in coordinate form are written 
as 



eo = A^ ^dt, ei = ei^S^ 



(2) 



In cosmological contexts it is customary to use the Hub- 
ble scalar H = 0/3 34]. From these it becomes obvious 
that in FRW geometries, i coincides with the scale fac- 
tor. Finally, it can be shown that these kinematical fields 
are related through ^, i34] 



CT^i^ := U{tiUy) + V(fj,u^) - -Qh^i, 



^flU 



i[^"l/] 



V[^u^]. 



(6) 

(7) 



Transiting to the synchronous temporal gauge, we can 
set N to any positive function of t, or simply A^ = 1. 
Thus, we extract the following restrictions on kinematical 
variables: 

< = diag(0,-2o-+,cr+,cr+), w^^ = 0, (8) 



The metric ^ can describe three geometric families, 
that is 

sin?? for k = +1, 

S{i}) = { 1} for fc = 0, (3) 

sinhi? for fc = — 1, 



where 



1 d 
3di 



In 



ei 

62 = 



(9) 



known respectively as Kantowski-Sachs, Bianchi I and 
Bianchi HI models. 



^ Covariant spacetime indices are denoted by letters from the sec- 
ond half of the Greek alphabet. 



Finally, note that the Hubble scalar can be expressed in 
terms of ei^ and e-i' as 



B 



1 d 



[lneii(e22) 



2n21 



(10) 



B. The action and dynamical variables 



Let us now construct /(i?) cosmology in such a geo- 
metrical background. In the metric formalism the action 
for /(i?)-gravity is given by [2lll2a| 



5„ 



dx" 



/(i?)-2A + £(") 



(11) 



where /(i?) is a function of the Ricci scalar i?, and £('"' 
accounts for the matter content of the universe. Addi- 
tionally, we use the metric signature (—1, 1, 1, 1), Greek 
indices run from to 3, and we impose the standard units 
in which c — 8ttG = 1. Finally, in the following, and 
without loss of generality, we set the usual cosmological 
constant A = 0. 

The fourth-order equations obtained by varying the 
action (fTTj) with respect to the metric write: 



^/iZ^ — 



T, 



(m) 



f'iR) 






(12) 



where the prime denotes differentiation with respect to 
R. In this expression T;)™ denotes the matter energy- 
momentum tensor, which is assumed to correspond to a 
perfect fluid with energy density p,„ and pressure Pm- 
Additionally, T/|, denotes a correction term describing a 
"curvature-fluid" energy-momentum tensor of geometric 
origin |35l. ISq: 



rriR _ 



1 



fiR) 



\g^, (/(i?) - RfiR)) 
+V^V,/'(i?)-5M,n/'(i?)], (13) 



where V^ is the covariant derivative associated to the 
Levi-Civita connection of the metric and D = VV^. 
Note that in the last two terms of the right hand side 
there appear fourth-order metric-derivatives, justifying 
the name "fourth order gravity" used for this class of 
theories ^2]. By taking the trace of equation ([121 and 
re-ordcring terms one obtains the "trace equation" (equa- 
tion (5) in section IIA of (sj) 



3n/'(i?) + RfiR) - 2/(i?) = T, 



(14) 



where T = T^ is the trace of the energy-momentum ten- 
sor of ordinary matter. 

In the phenomenological fluid description of a general 
matter source, the standard decomposition of the energy- 
momentum tensor T^i, with respect to a timelike vector 
field u^ is given by 



where /i denotes the energy density scalar, P is the 
isotropic pressure scalar, g^ is the energy current density 
vector (q^ u^ = 0) and tt^i^ is the trace-free anisotropic 
pressure tensor {tt^^u'^ = 0, tt^ = 0,7r^^ = tt^^). 

The matter fields need to be related through an ap- 
propriate thermodynamical equation of state in order to 
provide a coherent picture of the physics underlying the 
fluid spacetime scenario. Applying this covariant decom- 
position to the "curvature-fluid" energy-momentum ten- 
sor P^ we obtain 



fJ- 



P = — 



f(R)-Rf{R) + 6Hif'{R) 

r{R) 



-f{R) + RfiR) mj-J'iR) - 2i,f'{R) 



f'iR) 



(16) 



Finally, the anisotropic pressure tensor is given by tt^ 

diag(0, — 27r-|_,7r-(_,7r-|_). where 



f'iR) 



dt 



(17) 



C. The cosmological equations 

Let us now present the cosmological equations in 
the homogeneous and anisotropic Kantowski-Sachs met- 
ric. In particular, the Einstein's equations ()12|) in the 
Kantowski-Sachs metric read 

34 - 3i/2 -^K=-fi- -f" 



f'iRV 

f'iR) 

Pn 



- 3(a+ + HY - 2a+ - 2H - 'K = ^^ +P-2n+, 

P + TT+. 

(18) 



3al + 3a+H - 3H^ + n+ - 2H 



In these expressions pm and Pm are the energy density 
and pressure of the matter perfect fluid, and their ratio 
gives the matter equation-of-state parameter 



P,i 



(19) 



Furthermore, ^K is the Gauss curvature of the 3-spheres 
[3J] given by 



^K^{e^)^ 
and its evolution equation writes 



^k = -2{a+ + H) i^K) 



(20) 



(21) 



Additionally, the evolution equation for ei^ reads (see 
equation (42) in section 4.1 of [30l |) 



T„, 



p,u^t 



2q(f,u^ 



Ph„ 



(15) 



ei^ = -{H-2a+)ei^. 



(22) 



Finally, observing the form of the first of equations p8) 
one can define the various density parameters of the sce- 
nario at hand, namely [3M| the curvature one 



n. 



MP' 



the matter one 



^ ^tn — 



Pm 

the "curvature-fluid" one 

^ ^curv.fl 

and the shear one 



_p_ 



(23) 



(24) 



(25) 



(26) 



satisfying n^ + ^,n + ^curv.fi + ^i^ = 1. 

Now, the trace equation ([T^ in the Kantowski- Sachs 
metric reads 

-3-^fiR)-9Hj^f'{R)+Rf'{R)-2f{R) = -pm+ipm, 

(27) 
and the Ricci scalar writes 



R = 12H^ + 6al 



6H + 2^K. 



(28) 



At this stage we can reduce equations (fT8)) along with 
([27]) . with respect to H, a\ and "^K [S^, acquiring the 
Raychaudhuri equation 



H 



6f'{R) 



[Pm. + 3pm] — 



H d 
2f'{R)di 
3 d^ 



f'iR) 



fJR) 

,f'{R) .f'{R)dt' 



f'{R) 



(29) 



the shear evolution 



1 
6 



H' 



3/'(i?) 



R- 



m_ 

and the Gauss constraint 



^^1/'^-)^ (-) 



^K = 3al^3H' + ^^ + ^- 



R 



fiR) 



f'iR) 



3H d 
7W)di 



f'iR)- 

(31) 



It proves convenient to use the trace equation (P7)) to 
eliminate the derivative -^f'{R) in (1^^ . obtaining a sim- 
pler form of Raychaudhuri equation, namely 



H 



-H^ - 2al 



Pn 



fiR) 



' Hif'iR). 



3/'(i?) 6/'(i?) f'iR) dt 



(32) 



Furthermore, the Gauss constraint can alternatively be 
expressed as 



H 



1 



df'jR) 
dt 

2 f'iR) 



2/^ 2 , Pm 

3^ = "-^ + 3rM 



R- 



fiR) 



f'iR) 



df'(R) 
dt 

f'iR) 



(33) 



In summary, the cosmological equations of fiR)- 
gravity in the Kantowski-Sachs background are the "Ray- 
chaudhuri equation" (|32|) . the shear evolution (|30l) . the 
trace equation ([27| , the Gauss constraint ([33l) , the evolu- 
tion equation for the 2-curvature ^K (pij) . and the evolu- 
tion equation for ei^ (|22)) . Finally, these equations should 
be completed by considering the evolution equations for 
matter sources. 



D. i?"-gravity 

In this subsection we specify the above general cosmo- 
logical system. In particular, in order to continue we have 
to make an assumption concerning the function fiR)- We 
impose an ansatz of the form /(i?) = J?" [2I [22|. Isil. liot . 
since such an ansatz does not alter the characteristic 
length scale (and General Relativity is recovered when 
n — 1), and it leads to simple exact solutions which allow 
for comparison with observations [4l|, |42] . Additionally, 
following [33, [Saj we consider that the parameter n is re- 
lated to the matter equation-of-state parameter through 



n= -(l + w). 



(34) 



Such a choice is imposed by the requirement of the 
existence of an Einstein static universe in FRW back- 
grounds, which leads to a severe constraint in fiR), 
namely fiR) = 2A -h i?i(i+'") with w ^ -1. Inversely, 
it can be seen as the equation-of-state parameter of ordi- 
nary matter is a fixed function of rt, if we desire Einstein 
static solutions to exist. The existence of such a static 
solution (in practice as a saddle- unstable one) is of great 
importance in every cosmological theory, since it con- 
nects the Friedmann matter-dominated phase with the 
late-time accelerating phase, as required by observations 
[35I [36{ . Thus, if we relax the condition n = |(1 + w) 
in /J^-gravity, in general we cannot obtain the epoch- 
sequence of the universe. Finally, note that the constraint 
— 1/3 < ui < 1 that arises from the satisfaction of all the 
energy conditions for standard matter, imposes bounds 
on n, namely n e [1,3], with the most interesting cases 
being those of dust iw — 0, n — 3/2) and radiation fiuid 
(w = 1/3, n = 2). 

The cosmological equations for JS^-gravity, in the ho- 
mogeneous and anisotropic Kantowski-Sachs metric ([l}, 
are obtained by specializing equations ([32]) . ([30l) . ([33]) . 
([27| . (|2T|) . ([22I) and consider the conservation equation 
for standard matter. 



Equations (|32l), dSO)), dM]) and ^, become respec- 
tively 



H + H' 



-2crj 



3nR 



^"- + i^i?+(n-l)i/|,(35) 



i?' 



a+ = —(y^ — iHaj^ -)- H" 



inR" 



n — 1 



-i? 



R 



-{n-l)ia+-H)-, (36) 



i7 



n- li? 
2 ^ 



1, 



3 

n- 1 
6n 



ii = (Ji 



3ni?""i 

^2 n2 



-i? 



4 ;r2' 



(37) 



^ . ^-\R+'^^--\j^-3H^-in-2)^. 
R 3n(n-l) 3n(n-l)i?"-i i? ^ ' R^ 

(38) 
Additionally, we consider the evolution equation for the 
2-curvature ^K given by (PTj) . as well as the matter con- 
servation equation 



Pm = -2nHp,n- 



(39) 



Finally, in order to close the equation-system we consider 
also the propagation equation ([22l) . 



III. PHASE-SPACE ANALYSIS 

In the previous section we formulated /(i?)-gravity in 
the case of the homogeneous and anisotropic Kantowski- 
Sachs geometry, focusing on the f{R) = -R" ansatz. Hav- 
ing extracted the cosmological equations we can inves- 
tigate the possible cosmological behaviors and discuss 
the corresponding physical implications by performing 
a phase-space analysis. Such a procedure bypasses the 
complexity of the cosmological equations and provides 
us the understanding of the dynamics of these scenarios. 



A. The dynamical system 

In order to perform the phase-space and stability anal- 
ysis of the models at hand, we have to transform the 
cosmological equations into an autonomous dynamical 
system [43]. However, since the present system is more 



complicated, in order to avoid ambiguities related to the 
non-compactness at infinity we define compact variables 
that can describe both expanding and collapsing models 
pjll |44| . This will be achieved by introducing the auxil- 
iary variables [4^, |4^ : 



Q = 



H 



D 



(n-l)R {n-l)R 

x = — „„„ , y ^ . „^ , z 



K = 



2RD 



6nD^ 



3D^ 

where we have defined 



P 1 - "=1 



inR^-^D'^' 
(40) 



D = 



\ 



H 



n-lR 
2 R 



3 



'-^K. 



(41) 



Furthermore, we introduce the time variable r through 

D 



dr 



n-1 



dt. 



(42) 



and from now on primes will denote derivatives with re- 
spect to T. Finally, note that in the following we focus 
on the general case n 7^ 1. The dynamical investigation 
of Kantowsky-Sachs model in General Relativity (that is 
for n — 1) has been performed in [47]. 

In terms of these auxiliary variables the Gauss con- 
straint ([57)1 becomes 



x^ + y + z + T.^ = 1. 
Moreover, the D-definition (|^T|) becomes 

{Q + xf + K = l. 



(43) 



(44) 



Finally, from the definitions (|40l) we obtain the bounds 
y > 0, z > and K > 0. Therefore, we conclude that 
the auxiliary variables must be compact and lie at the 
intervals 



Qe 1-2,2], 
xe 1-1,1], 
Ke[0,i], 



Se 1-1,1], 

ye 10,1], zelO,i] 



(45) 



while Ei^ is un-constrained. 

In summary, using the dimensionless auxiliary vari- 
ables (I40p , along with the two constraints (H51) and (|44l) , 
we reduce the complete cosmological system to a five- 
dimensional one given by: 



/ = (1 - n)j:Q^ -{n-1) [-3x^ + 2S.t + {n - 3)2^ + n {x^ + y - l) + l] Q^ 
-{n - 1) {x [Sx^ + (x - 3S)S + n (a;2 + I]2 + y - l) - l] - E} g 
~x^ + E^ + n (a;^ - I]2 + y - 1) + 1, 



(46) 



E' = -(n - 3)(n - 1)(Q + x)T,^ - {n - 1){Q + x- 1)(Q + x + 1)Y.^ 

-{n - 1)(Q + x) [{n - 3) {x^ - l) + ny] E + (n - 1) [(Q + xf - l] 



(47) 



x' = -(n - 3)(n - l)a;^ - (ri - 1) [(« - 3)0 + E] x^ - (n - 1) [-3^2 + 2gS + n (E^ + y - 2) + 5] a;^ 

-(n - 1) {Q [{Q ~ 3S)S + n (S^ + y - l) + 3] - S} ,t + (n - 2) [-S^ + n (S^ + y - l) + l] , (48) 

y' = -2(n - l)n{Q + x)y^ + y {-2 [{n ~ 4)n + 3] (Q + .x)E2 - 2(n - l)(g + a; - l)(g + a; + 1)1] 

-2 {(g + x) (a;2 - 1) n2 + (-4x3 _ 4q^2 + 2a; + Q) n + a;[3a;(g + a;) - 2]}} , (49) 



(^i^)' = Ex^{n ~ 1) {~[n - 3){Q + x)J:^ - [{Q + xf - 3] S + (-g - x) [-3x^ + n (x'^ + y - l) + l]} . (50) 



r 



Proceeding forward, as we observe, equation (|50p de- 
couples from (|46l) - (l49l) . Thus, if one is only interested in 
the dynamics of Kantowsky-Sachs /(i?)-models, all the 



J 



relevant information of the system is given by the remain- 
ing equations. Therefore, we can reduce our analysis to 
the system (|l^ -(|49 p . defined in the compact phase-pace 



* = {a 



S' < l,|Q + a;| < l,a;G [-1,1], yG [0,1], EG [-1,1], ge [-2,2]}. 



(51) 



r 



Finally, while for the purpose of the present work it 
is adequate to investigate the system PS)) - (P^ . leaving 
outside the decoupled equation ([50]) , for completeness we 
study equation (|50| in Appendix |B] 



B. Invariant sets and critical points 

Using the auxiliary variables (PO)) the cosmological 
equations of motion p5p - ([5^ . were transformed into the 
autonomous form (|^5|) - ([5(I)) . that is in a form X' — f(X), 
where X is the column vector constituted by the auxil- 
iary variables and f(X) the corresponding column vector 
of the autonomous equations. The critical points Xc are 
extracted satisfying X' = 0, and in order to determine 
the stability properties of these critical points we expand 
around Xc, setting X = Xc-I-U with U the perturbations 
of the variables considered as a column vector. Thus, up 
to first order we acquire U' = H • U, where the matrix 
H contains the coefficients of the perturbation equations. 
Therefore, for each critical point, the eigenvalues of H de- 
termine its type and stability. In particular, eigenvalues 
with negative (positive) real parts correspond to a stable 



(unstable) point, while eigenvalues with real parts of dif- 
ferent sign correspond to a saddle point. Lastly, when at 
least one eigenvalue has zero real part, the corresponding 
point is a non-hyperbolic one. 

Now, there are several invariant sets, that is areas of 
the phase-space that evolves to themselves under the dy- 
namics, for the dynamical system P5|) - (H^ . There are 
two invariant sets given by g-t-a; = ±1, corresponding to 
K = 0, that is 62^ /D = 0. However, the solutions of the 
evolution equations inside this invariant set do not corre- 
spond to exact solutions of the field equations, since the 
frame variables has to satisfy det(eQ) 7^ 0. Nevertheless, 
it appears that this invariant set plays a fundamental 
role in describing the asymptotic behavior of cosmologi- 
cal models. 

Another invariant set is y = 0, that is pm = 0, i? = 
12H^ + 6H+6(7^ + 2'^K = 0, provided D is finite, or p,„ = 
0, R = Rq, with Rq a constant, provided D — > 00. This 
invariant set contains vacuum (Minkowsky) and static 
models. 

Another invariant set appears in the case of radiation 
background (n = 2,w — 1/3), namely that of x = 0. 

Finally, we have identified the invariant set x"^ + y + 



S^ = 1, which contains cosmological solutions without 
standard matter. 



C. Local analysis of the dynamical system 

The system (f46 |) - (|49|) admits two circles of critical 
points, given by Cc : E^ + Q^ — 2e(3 — 0, x + Q — e, y — 
(we use the notation e — ±1). They exist for n E [1,3], 
they are located in the boundary of ^, and they cor- 
respond to solutions in the full phase space, satisfying 
K = y = z = 0. In order to be more transparent, let us 
consider the parametrization 



Q = e + sin u 
Cf -.^ { S = cosu , u e [0, 27r] 
X = — sinw. 



(52) 



For all the critical points we calculate the eigenvalues 
of the perturbation matrix H, which will determine their 
stability. These eigenvalues, evaluated at C^, are 



{0, -2(n-l)(cosu-2e), 
-2(n - 1) [{n - 2) sinu - e(3 - n)] 
2(n-2)sinii + 6e(n- 1)}. 



(53) 



For the values of u such that there exists only one 
zero eigenvalue, the curves are actually "normally hyper- 
bolic" ^ , and thus we can analyze the stability by analyz- 
ing the sign of the real parts of the non-null eigenvalues 
p8| . Therefore, we deduce that: 



our notation: £_ = C_|„=3^/2, A/V = C+|„=3^/2 and 
A/"- = C_|„=^/2,£+ = C+L ,r/9 - Moreover, and con- 
trary to the investigation of [36|, we obtain four new 
critical points, which are a pure result of the anisotropy. 
They are labelled by P^ := (Q = e, E = e, x = 0, y = 0) 
(Q = e, S = — e, .T = 0, J/ = 0), and obviously 



and P^ 

C-\ 



C+\u=OtPi — C-\u=TT and Pj 
Thus, it is easy to see that: 



C. 



P~ 

-TT; -* 2 



• For — i < w < I, C+ is unstable and £_ is sta- 
ble. The critical point N+ (A/-) is always unstable 
(stable) except in the case w = —^. These results 
match with the results obtained in |36| . 

• The critical points P^ and P2 exist always. They 
are non- hyperbolic, presenting an 1-dimensional 
center manifold tangent to the line x + Q = 0. 
P^2 (^1^2) have a 3D stable (unstable) manifold 



provi' 



ded 



'\ <w <l. 



The curves of critical points Ce, and the representative 
critical points A^, £e, Pf and P2^ are enumerated in 
Table D 

Until now we have extracted and analyzed the stability 
of the curves of critical points Ce of the system (|151) - (|15)) . 
and their representative critical points. However, the 
system (|46p - ([49|) admits additionally ten isolated criti- 
cal points enumerated in Table |lTl where we also present 
the necessary conditions for their existence. ^ 

The eigenvalues of the Jacobian matrix are presented 
in Appendixl^and thus here we straightway provide each 
point's type. 



• No part of C+ (C_) is a stable (unstable). Thus, 
any part of C+ (C_) is either unstable (stable) or 
saddle. 



For 



< w < 



< 



^ ^ ^ ^ i one branch of C_ ( ), '^\ ' 

3 — 6' ^ 3u? — 1 

sin(u) < 1) is stable, and one branch of C+ 
^3(3«^+il ^ gij^(_y) < i^ jg unstable. 



3to-l 



• For 



< w < 



g ^ ^ ^ 3 

whole C+ is unstable 



the whole C_ is stable and the 



• For I 



< w < 1, one branch of C_ ( 



3(w-l) 
3t«-l 



< 



sin(u) < 1) is stable, and one branch of C+ 
( 3j«^i "^ sin(— u) < 1) is unstable. 

Note that amongst the curves Ce, we can se- 
lect the representative critical points labelled by 
^f^ := (g = 0, E = 0, a; = e, y = 0) and C^ ■= 
(g = 2e, E = 0, X = -e, y = 0) described in [11]. In 



The two critical points A^ exist for —h<w< 
1. The critical point A+ {A-) is a sink, that is a 
stable one (source, that is an unstable) provided 
■u;+ < w < 1, where w+ = ^ (-2 -h Vs) « -0.09. 
Otherwise they are saddle points. Equivalently, we 
can express the stability intervals in terms of n, 
using the relation (l34t . Thus, the aforementioned 
interval becomes n-|_ < n < 3, where n+ = 3(w+ + 
l)/2« 1.37. 



• The critical points Be exist for — i < w < |. 
Thus, they are always saddle points having a 2- 
dimensional stable manifold and a 2-dimensional 
unstable manifold if 



i<w<l 



• The critical points P| exist for — i < w < w^^ . The 
are non-hyperbolic, and due to the 1-dimensional 
center manifold presented in the Appendix \^ the 
stable (unstable) manifold of P^ {P3) is 3D for 
— i < w < w+. 



^ Since we are dealing with curves of critical points, every such 
point has necessarily at least one eigenvalue with zero real part. 
"Normally hyperbolic" means that the only eigenvalues with zero 
real parts are those whose corresponding eigenvectors are tangent 
to the curve |48||. 



^ Strictly speaking, the system I I46I I- II49I I admits two more critical 
points with coordinates Q = ±2, S = ±1, x = 0,y = 0. However, 
they are unphysical since they satisfy | Q + a:: | > 1 and thus ^K = 



Cr./Curve 


Q 


E 


X 


y Existence 


Stability 


A/; 








1 


always 


unstable 


M- 








-1 


always 


stable 


C+ 


2 





-1 


always 


unstable for 1.25 < n < 2.5 


C- 


-2 





1 


always 


stable for 1.25 < ti < 2.5 


Pt 


1 


1 





1 < n < 3 


non-hyperbolic with 3Z) unstable manifold 


Pi 


-1 


-1 





l<n<3 


non-hyperbolic with 3D stable manifold 


P^ 


1 


-1 





1 < n < 3 


non-hyperbolic with 3D unstable manifold 


P2 


-1 


1 





l<n<3 


non-hyperbolic with 3D stable manifold 


c+ 


1 + sin u 


COS It 


— sin It 


always 


unstable for 1.25 < n < 2.5 


c_ 


— 1 + sin u 


COSM 


— sin It 


always 


stable for 1.25 < n < 2.5 



TABLE I: The curves of critical points Cj, and the representative critical points Mt, £e, Pi and PI of the cosmological system 
(we use e = ±1). u varies in [0, 27r]. For more details on the stability of Ce see the text. 



Cr. P. 


Q 


E 


X 




y 


Existence 




Stability 


A+ 


2n-l 
3(n-l) 





n-2 
3(n-l) 




(2n-l)(4n-5) 
9(n-l)^ 


1.25 < n< 3 




stable for rt+ < rt < 3 
saddle for 1.25 < n < n+ 


A- 


2n-l 
3{n-l) 





n-2 
3{n-l) 




(2n-l)(4n-5) 
9(n-l)^ 


1.25 < n< 3 




unstable for n+ < n < 3 
saddle for 1.25 < n < n+ 


13+ 


i 
3-n 





n-2 
n-3 







l<n< 2.5 




saddle 


B- 


i 
3-n 





n-2 
n-3 







1 < n< 2.5 




saddle 


P+ 


i 
2-n 





n-i 
2-n 




n-i 


1 < n < n+ 


non 


-hyperbolic with 3D stable manifold 


n(n-2)^ 


P3- 


i 
2-n 





n-i 
2-n 




n— i 


1 < n < n+ 


non 


-hyperbolic with 3D unstable manifold 


n(n-2)^ 


Pt 


2n2-5n + 5 
7n^-16n + 10 


2n2-2n-l 
7n^-16n + 10 


3(n-l)(n-2) 
7n^-16n + 10 


9(4: 

1 


n,"'-10n-|-7)(n-i)"' 
;7n'i-16n + i0)-^ 


n+ < n < 3 




saddle with 3D stable manifold 


Pa 


2n^-5n-|-5 
7n^-16n + 10 


2n^-2n-l 
7n^-16n-|-10 


3(n-l)(n-2) 
7n^-16n-f-10 


9(4: 

1 


n,^-10n + 7)(n-i)^ 
;7n^-16n + 10)^ 


n+ < n < 3 




saddle with 3D unstable manifold 


Pt 


1 
n — 2 


V2n-5 
n-2 


n-3 
n-2 







2.5 <n < 3 


non 


-hyperbolic, with a 2D center manifold 


P^ 


1 
n-2 


V2n-5 
n-2 


n-3 
n-2 







2.5 < n < 3 


non 


-hyperbolic, with a 2D center manifold 



TABLE IL The isolated critical points of the cosmological system. We use the notation n+ 



1+73 



1.37. 



• The critical points PI exist for w+ < w < 1. P^ 
{P4) has a 3D stable (unstable) manifold and a ID 
unstable (stable) manifold provided w+ < w < 1. 
Thus, they are always saddle points. 

• The critical points P| exist for | < w < 1. They are 
non- hyperbolic, there exists a 2D center manifold 
and Pc^ {P5) has a 2D unstable (stable) manifold. 

Lastly, as we have mentioned, in this subsection we 
have investigated the system (|46|) - (|49|) . leaving outside 
the decoupled equation (|50p . However, for completeness, 
we study equation (|50p , and especially the stability across 
the direction Ei^, in Appendix [Bl 



D. Physical description of the solutions and 
connection with observables 



Let us now present the formalism of obtaining the 
physical description of a critical point, and also connect- 
ing with the basic observables relevant for a physical dis- 
cussion. These will allow us to describe the cosmological 
behavior of each critical point, in the next section. 



Firstly, around a critical point we obtain first-order ex- 
pansions for ei^, ^K, pm and R in terms of r, considering 
the versions of equations ((22)) . ((2T|) and ((39)) . respectively 
given by 



(ei^)' 


= -(7i-l)[Q*-2S*]ei\ 


W 


= -2(n-l)[0*+I]*]X 


P',n 


= -2n(rt-l)gV™, 


R! 


= 2x*R, 



(54) 



where the star-upperscript denotes the evaluation at a 
specific critical point, and the prime denotes derivative 
with respect to r. The last equation follows from the 
definition of x given by (j40p . In general we can consider 
the case x* ^ G and y* ^ 0, since in the simple case 
of X* =0 and y* 7^ we obtain R = Rq, with Rq a 
constant, while in the case y* = we acquire R = 0. 

In order to express the above determined functions of 
T in terms of the comoving time variable t, we invert the 
solution of 



dt 
1^ 



1 



D* 



(55) 



with D* being the first-order solution of 




D' = D{n - 1)T^, 




where 




T* ^x* -j:* + iQ*+x*) 


-3a;*Vs*x* 




+ (Q*-3E*)S]*+n( 


x*^ + X*^ + y* - 


-01 



(56) 



.(57) 



Solving equations (j55p - (|56p (with initial conditions 
D{0) = Da and t{0) = to) we obtain 



t{r) 



,(l-n)rT* 



^oT* 



to. 



(58) 



Thus, inverting the last equation for r and substituting 
in the solution of (l54l) with initial conditions ei^(O) = 
eij, ^^(0) = ^Kq, p.,„(0) == Pmo,R{0) = i?o, we acquire 



ei\t) ^ eil{Do{to ~ t)T* + 1) 

Prnit) ^ PmoiDoito - t)T* + 1) 



Q' -2S' 



2(Q*+S'') 



2nQ' 



R{t) = i?o(i^o(io - t)T* + 1)0^^ 



(59) 



Finally, it can be shown that the length scale i along the 
flow lines, defined in ([5]), can be expressed as [34] 



Finally, the various density parameters defined in ((23)) - 
(|26p , in terms of the auxiliary variables straightforwardly 
read 



n,.. 



{Q + xf -1 



^ ^curv.fl 



l-x"^ -y-Y? 
y - 2xQ 



Q 



(65) 



In conclusion, at each critical point we can calculate 
the values of the basic observables q, WeS and the various 
density parameters, and also calculate the specific phys- 
ical solution, that is obtain the behavior of f(t), pm{t) 
and R{t). 

Finally, it is interesting to notice that in Kantowski- 
Sachs geometry one can easily handle isotropization. In 
particular, the geometry becomes isotropic if ct+ becomes 
zero, as can be seen from ([1]) and (|9]). Thus, critical 
points with E = (or more physically fla- = 0) corre- 
spond to Friedmann points, that is to isotropic universes, 
and when such an isotropic point is an attractor then we 
obtain asymptotic isotropisation in the future [3l| . 



e^£oiii-tT*)- 



(60) 



where Iq = [(dj) {^Kq)] ^ and h = Dq^T* -f 1. In 
summary, expressions (|59p and (j60p determine the so- 
lution, that is the evolution of various quantities, at a 
critical point. 

Let us now come to the observables. Using the above 
expressions, we can calculate the deceleration parameter 
q defined as usual as [3J| 






(61) 



Additionally, we can calculate the effective (total) 
equation-of-state parameter of the universe WcS, defined 
conventionally as 



WcS- 



Ptot 
Ptot 



f'iR) 



+ p 



f'{R) 



M 



(62) 



where Ptot and ^tot are respectively the total isotropic 
pressure and the total energy density as they can be read 
from equation ([TSl) . where P and p, are given by ([T6| . 
Therefore, in terms of the auxiliary variables we have 



x{2Q + x) + l 



ny 



(n-l)02 



(63) 



WeS 



'2ny 1 , ^ 
1 ^ _ (64) 

3(n-l)(a;2 + 2Qa; + I]2-l) ^3 ^ ' 



IV. COSMOLOGICAL IMPLICATIONS 

In the previous sections we formulated /(-R)-gravity in 
the case of the homogeneous and anisotropic Kantowski- 
Sachs geometry, focusing on the f{R) — i?" ansatz, and 
we performed a detailed phase-space analysis. Thus, in 
the present section we discuss the physical implications 
of the mathematical results, focusing on the physical be- 
havior and on observable quantities. 

First of all, for each critical point of Tables HI and HIl we 
calculate the effective (total) equation-of-state parame- 
ter of the universe Wcs using (|64|) . and the deceleration 
parameter q using (p5)) . and the values of the density 
parameters flk, ^m, ^curv.fi, and fta using ([S5)) . These 
results are presented in the Table IIIII 

Furthermore, for each critical point we use (I59|) and 
([60)) . in order to extract the behavior of the physically 
important quantities £(t), pm{t) and R{t) at this criti- 
cal point. As we have mentioned (,{t) is the length scale 
along the flow lines and in the case of zero anisotropy (for 
instance in FRW cosmology) it is just the usual scale fac- 
tor. Additionally, Pm{t) is the matter energy density and 
R{t) is the Ricci scalar. These solutions are presented in 
the last column of Tables IIVI and [Vj 

Finally, in the last column of Tables IIVI and |V] we also 
present the physical description of the corresponding so- 
lution, taking into account all the above information. In 
particular, since the auxiliary variable Q defined in (|40p 
is the Hubble scalar divided by a positive constant, Q > 
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Cr.P. 


Weft- 


<? 


fife 


f^m 


f^curv.fl 


Q.a 


A± 


6n''-7n-l 


2n^-2n-l 
(n-l)(2n-l) 








1 





3(n-l)(2n-l) 


B± 


3 


1 





5-2n 


2(n - 2) 





P^.P^ 


1 
3 


2 











1 


Pt 


1 








-2n + 2 + ^ 


2n-l-i 





Pt 


n(4(9-4n)n-21)-l 
3n(8(n-3)n+21)-3 


6-3n 1 


3(2n^ -2n-l) (4n^ -lOn + 7) 





3(n-l)(8n^-24n2+21n-l) 
(2n2-5n + 5)^ 


(2n2-2n-l)^ 


n(2n-5)+5 ^ 


Pt 


i 
3 


2(n-2) 








2(3 - n) 


2n-5 


c± 


1 
3 


2 








2 sin{u) 
sin{u)±l 


cos^(u) 


l±sin(u) 


(siii(u)±l)^ 



TABLE III: The values of the basic observables, namely the effective (total) equation-of-state parameter Weff , the deceleration 
parameter 5, the curvature density parameter n^, the matter density parameter f2m, the "curvature-fluid" density parameter 
ficurv.fl and the shear density parameter fio-, defined in (I63p -(|65 p . at the isolated critical points and curves of critical points 
of the cosmological system. We display the information about the non-isolated critical points P^^P^ to emphasize they are 
solutions dominated by shear. 



corresponds to an expanding universe, while Q < to a 
contracting one. Furthermore, as usual, for an expanding 
universe g < corresponds to accelerating expansion and 
g > to decelerating expansion, while for a contracting 
universe g < corresponds to decelerating contraction 
and g > to accelerating contraction. Additionally, if 
Weflf < — 1 then the total equation-of-state parameter of 
the universe exhibits phantom behavior. Lastly, critical 
points with E = correspond to isotropic universes. 

Let us analyze the physical behavior in more details. 
The critical point Aj^ is stable for 1.37 < n < 3, and 
thus at late times the universe can result in it. It corre- 
sponds to an isotropic expanding universe, in which the 
expansion is accelerating. Additionally, since the curva- 
ture energy density is zero, this point corresponds to an 
asymptotically flat universe (this is equivalent to close 
or open FRW universes which may possess zero curva- 
ture solutions as asymptotic states). Furthermore, in 
the case 2 < n < 3 the total equation-of-state param- 
eter of the universe exhibits phantom behavior. Thus, 
the anisotropic Kantowsky-Sachs _R"-gravity can lead the 
universe not only to be accelerating, but also to be in the 
phantom "phase" , a result of great cosmological interest. 
We have to mention here that since point A+ is asymp- 
totically an FRW one, one expects to obtain phantom 
behavior in isotropic (FRW) _R"-gravity too. Although 
in the phase-space analysis of FRW i?"-gravity [Sg] , the 
authors have focused only on acceleration, without ex- 
amining the total equation-of-state parameter of the uni- 
verse, it is easy to see that if ones defines it and exam- 
ines its features he will find phantom behavior in that 
case too. This is also the case in Bianchi I and Bianchi 
III _R"-gravity [3l|, |43], where the authors would have 
found the phantom behavior if they had calculated the to- 
tal equation-of-state parameter. Therefore, we conclude 
that the phantom behavior is a result of the modification 
of gravity, as it has been discussed in detail in the lit- 
erature |2l| - |23j . Finally, in the case where the matter is 
dust {w = 0, that is n = 3/2), Pmit) behaves like ^(i)~^ 
as expected, and this acts as a self-consistency test for 
our analysis. Additionally, note that in the special case 



n = 2, that is when w = 1/3, i.e when radiation dom- 
inates the universe, the aforementioned stable solution 
corresponds to a de-Sitter expansion (this is not a new 
feature since de-Sitter solutions are known to exist in 
Bianchi I and Bianchi III i?"-gravity [3l|). This is of 
great significance since such a behavior can describe the 
inflationary epoch of the universe. 

In the above critical point the isotropization has been 
achieved. Such late-time isotropic solutions, than can at- 
tract an initially anisotropic universe, are of significant 
cosmological interest and have been obtained and dis- 
cussed in the literature [ij, |31| . The acquisition of such 
a solution was one of the motives of the present work, 
as of many works on anisotropic cosmologies, since, as 
we discussed in the Introduction, it is the only robust 
approach in confronting isotropy of standard cosmology. 
The fact that this solution is accompanied by acceleration 
or phantom behavior, makes it a very good candidate for 
the description of the observable universe. 

The critical point A- which corresponds to an 
isotropic contracting universe is not stable and thus it 
cannot be the late-time state of the universe. Similarly, 
the points B+ and B- , which correspond to isotropic ex- 
panding and contracting universes respectively, and in 
which the total matter/energy mimics radiation, are sad- 
dle points and thus they cannot be the late-time solution 
for the universe, too. The points P^, which correspond 
to decelerating expansion (for e = +1) and to accelerat- 
ing contraction (for e = —1), are non- hyperbolic with a 
2D center manifold, and thus, generically, the universe 
cannot be led to them. 

The critical points P^' and P2~ which correspond to ac- 
celerating contraction, in which the total matter/energy 
behaves like radiation, are non-hyperbolic critical points, 
but they do possess a 31? stable manifold. By an explicit 
and straightforward computation of their center mani- 
folds [50] we deduce that for 2 < n < 3 each center man- 
ifold is locally asymptotically stable (the equation gov- 
erning the dynamics on the center manifold is a gradient- 
type differential equation with potential function having 
a degenerate local minimum at the origin), and thus Pj^ 
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Cr.P. T* Solution/description 



-nz 



n 



/ {1.25, 2} 



3(n-l)^ 



-4+ ^t^ m = I 'f'njZT lt\, PAt) = P^o [i] "" , R{t) = '"tr^ n = 2 

'- " ^ ^ [ n = 1.25 

Isotropic. Expanding. Accelerating for n+ < n < 3. Phantom for 2 < n < 3. 

For n = 2 exhibits dS behavior 



(fn(f,-tT*Y^ 77^2 r 1 -2n f (£i-tT*)^ h/ {1.25,2} 



n = 1.25 

Isotropic. Contracting. Decelerating for n+ < n < 3. Phantom for 2 < n < 3. 

For n = 2 exhibits collapsing AdS behavior 
Isotropic. Expanding. Decelerating. Total matter/energy mimics radiation. 

m = f„\/(^i - iT*), p„(f) = p„o(^i - tT*)-", i?(t) = 

Isotropic. Contracting. Accelerating. Total matter/energy mimics radiation. 



2 
n — 3 



"R;^ ^i5 ^(i) = ^o(^i - tT*) = ai + a2i, p™(i) = p™o(^i - tT*)""", R{t) = (^^4%=; 

Isotropic. Expanding. Zero acceleration. 

Isotropic. Contracting. Zero acceleration. 

p+ 5^-2) . . _ r ^o(^i -tT')^^ n / 2 ^ (,._^ li]-'" r.(,. _ f (I7^i%^ n / 2 



Expanding. Accelerating for n+ < n < 3. Phantom for 2 < n < M+. 
For n = 2 exhibits dS behavior 



^- 3(n-2) Z7. (£o(£i-tr*r^ n/2 ^ r,i-2" p.,. r (. _?^.). nfy 



Contracting. Decelerating for n+ < n < 3. Phantom for 2 < n < M+. 
For n = 2 exhibits collapsing AdS behavior 



n+ 


(3-2u) 
n-2 


Expanding. Decelerating. Total matter/energy mimics radiation. 


n^ 


(3-2n) 
n-2 


Contracting. Accelerating. Total matter/energy mimics radiation. 



TABLE IV: The behavior of £{t) (length scale along the flow lines), of pmit) (matter energy density) and of R{t) (Ricci scalar) 

l)(2n-l) _ 2n^-5n + 5 

n-2 ' *2 3(n-2) ' 



at the critical points of the cosmological system. We use the notations Si = —— — ^^^y — -, S2 = — ^"„, ^o\^ i ^s — y-^ 



P^(«) = &S^)^ n+ = i^ « 1.37 and M+ = i (5 + V2T) ^ 2.40. 



Cr.P. T* Solution/description 

'^ ^ r 1 —'Zn 

C+ -[3 + sin(lt)] m=£o{£l-ti:'Y+^''\ prr.it) =Prr.o[j^\ , R{t) = 

Expanding. Decelerating. Total matter/energy mimics radiation. 



C- 



-[-3 + sin(u)] £(i) = £o(^i-iT*)''-("',p„(i) = p„o[^] " , R{t) = 

Contracting. Accelerating. Total matter/energy mimics radiation. 



TABLE V: Physical behavior of the solutions at the curves of critical points of the cosmological system. We use the notation 

6+sin(n) 
3€+sin(n) 



/ \ 6+sinfn) 

Pe{u) - -^ ^—i- 
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and P2 are locally asymptotically stable [50|, and can 
attract the universe at late times. On the contrary, for 
I < n < 2, P^ and P2 are locally asymptotically un- 
stable (of saddle type). In summary, although P^ and 
P2 are not stable, they do have a significant probability 
to be a late-time state for the universe (this is realized 
for initial conditions on its stable manifold), or at least 
the universe can stay near these solutions for a long time 
before the dynamics remove it from them. On the other 
hand, the points P^ and P2 which correspond to decel- 
erating expansion, possess a 3D unstable manifold, and 
thus they cannot be a late-time solution of the universe. 

Non-hyperbolic critical point with a ZD stable mani- 
fold is also the critical point P^ , which corresponds to 
an asymptotically fiat isotropic expansion with zero ac- 
celeration, and thus it has also a large probability to be 
a late-time state of the universe. However, note that this 
solution corresponds to zero acceleration, and thus (.{t) 
is a linear function of t. On the other hand, the critical 
point Pg" (isotropic, contracting with zero acceleration) 
possesses a 3D unstable manifold and thus it cannot at- 
tract the universe at late times. 

The critical point P^ is saddle with a 3D stable mani- 
fold, and thus it has a large probability to be a late-time 
state of the universe. It corresponds to a non-flat, accel- 
erating expansion, for 1.37 ^n < Z, and furthermore in 
the case 2 < n < 3 it exhibits phantom behavior. Finally, 
for n ~ 2, that is for w = 1/3, it corresponds to a de- 
Sitter expansion. On the other hand P^ (decelerating 
contraction) is highly unstable and therefore it cannot 
attract the universe at late times. 

We mention here that the critical points Pf ,P2 ,P§ ,Pl 
and P| are not present in isotropic (FRW) i?"-gravity, 
as compared with [Sg . They arise as a pure result of the 
anisotropy, and this shows that the much more compli- 
cated structure of anisotropic geometries leads to radi- 
cally different cosmological behaviors comparing to the 
simple isotropic scenarios. 

Additionally, we have to analyze the behavior of the 
curves of critical points C^ . The points C_ , which cor- 
respond to accelerating contraction, are stable if 1.25 < 
n < 2.5, and thus they can be late-time solutions of the 
universe. On the other hand, C+, which correspond to 
decelerating expansion, are unstable and thus they can- 
not attract the universe at late times. 

Let us finish the physical discussion by referring to 
static solutions, in order to compare to the FRW case of 
|36| . From the cosmological point of view, static solu- 
tions possess £{t) = const., that is £{t) = and i{t)=0, 
in order to obtain H{t) = and H{t) = (note that one 
needs both conditions, since in a cosmological bounce or 
turnaround, that is when a universe changes from con- 
tracting to expanding or vice versa, £{t) is zero instantly, 
but then it becomes positive or negative again). Thus, 
using our auxiliary variable Q defined in (I40p static solu- 
tions should have Q = and Q = 0. In conclusion, since 
in all the aforementioned critical points Q = by defini- 
tion, static solutions are just those with Q — 0. At this 



point there is another important difference comparing to 
the isotropic case of |36|]. In particular, in FRW geom- 
etry, the Ricci scalar is R = 6{H + 2H^), and therefore 
static solutions have R = (however the inverse is not 
true, that is not all solutions with R = correspond to 
static ones, since one can have R — but with H and H 
non-zero). On the other hand, in the anisotropic case R 
is given by dH]), i.e P = 12ij2 + Qal +6H + 2'^K, that is 
we have also the presence of additional terms, and there- 
fore static solutions do not correspond to i? = unless 
CT+ and ^K are also zero, which is not fulfilled in general. 
Thus, in the anisotropic case one cannot use R, or equiv- 
alently the auxiliary variable y defined in (1401) . in order 
to straightway determine the static solutions, in contrast 
to the isotropic case |36| where the authors use y = for 
such a determination. 

Now, in order to present the aforementioned results in 
a more transparent way, we perform a numerical elabo- 
ration of our cosmological system, using a seventh-eighth 
order continuous Runge-Kutta method with absolute er- 
ror 10~^, and relative error 10~* [5l[. In Fig. [1] we 
depict some orbits in the invariant set y = 0, in the 
case of dust matter {w ^ 0,n ~ 3/2). As we observe, 




FIG. 1: Projection of the phase space on the invariant set 
y — 0, in the case of dust matter (w = 0, n = 3/2). The 
critical points Qf 2 have coordinates Qf :— (Q = 4e/3, E = 
2V2e/3,-e/3) and Q| — (Q = e/3, E = v^Se/S, 2e/3). All 
the points in the circle C+ are unstable whereas all the points 
in the circle C_ are stable. The heteroclinic sequences reveal 
the possible transition from expansion to contraction and vice 
versa (see text). 



we have the appearance of four critical points with co- 
ordinates Q\ := (Q = 4e/3,S = 2V2e/3, -e/3) and 
Q2 •= {Q = e/3, S = -\/5e/3, 2e/3), located in the in- 
variant curves C^. All the points in the circle C+ are 
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unstable whereas all the points in the circle C_ are sta- 
ble. Note that the critical points P| coincide with AC- 
In the figure we also display heteroclinic sequences [3M| 
of types 




(66) 

Thus, we can have evolutions in which the + branch 
and — branch are connected, that is we can have the 
transition from expansion to contraction and vice versa. 
This is just a cosmological turnaround and a cosmologi- 
cal bounce, and their realization in the present scenario 
reveals the capabilities of the model. It is interesting to 
note that such behaviors can be realized in FRW i?"- 
gravity [40, [S^] , however they are impossible in General 
Relativity Kantowsky-Sachs cosmology [45|. Therefore, 
we conclude that they are a result of the i?" gravitational 
sector and not of the anisotropy. 

In Fig. [2]we display some orbits in the case of radiation 
(w = 1/3, n = 2), where as we mentioned in subsection 
IIIIBI the invariant set x = appears. Particularly, we 
observe the existence of heteroclinic sequences of types 



A- 

B+ 
Pi 
P2 
P^ 



Pi 



p: 



A, 



A+ 
^Pt 
>Pi 



P^ 
P^ 



P. 



1 ' 



(67) 



revealing the realization of a cosmological bounce or a 
cosmological turnaround. Similarly to the isotropic case 
(see figure 5 of [30|) there is one orbit of type B+ — s- A+ 
and one of type A- — >■ B- . However, in the present case 
we have the additional existence of an heteroclinic se- 
quence of type A- — > P4 — > P4 — > A+ , correspond- 
ing to the transition from collapsing AdS to expanding dS 
phase, that is we obtain a cosmological bounce followed 
by a de Sitter expansion, which could describe the infla- 
tionary stage. This significant behavior is a pure result 
of the anisotropy and reveals the capabilities of the sce- 
nario. Lastly, note the very interesting possibility, of the 
eternal transition P^ — > P^ — > P^ — > P^ ■ ■ ■ which 
is just the realization of cyclic cosmology [5j. Bouncing 
solutions are found to exist both in FRW i?"-gravity [5JI, 
as well as in the Bianchi I and Bianchi III i?"-gravity [31| 
(see also [53|), and thus they arise from the i?" gravita- 
tional sector. However, in the present Kantowsky-Sachs 




FIG. 2: Projection of the phase space on the invariant set 
2; = 0, in the case of radiation {w — 1/3, n — 2). There is one 
orbit of type B+ — > .4+ and one of type A~ — > S_ . The 
existence of an heteroclinic sequence of type A- — > P^ — > 
P^ — > A+ corresponds to a cosmological bounce. 



geometry cyclicity seems to be realized relatively easily 
(however without accompanied by isotropization), that 
is without fine-tuning the model-parameters, which is an 
advantage of the scenario. 

Finally, in Fig. |3] we have extracted orbits located at 
a; = J/ = 0, in the case of radiation. The shaded region 
corresponds to the unphysical portion of the phase plane 
(note however that this region is not invariant, since an 
open set of orbits enter/abandon the unphysical bound- 
ary, and thus such evolutions have to be excluded too). 
As we observe, in this figure the last two heteroclinic 
sequences of ([S7)) are displayed. 

Let us make a comment here by referring to the cosmo- 
logical epoch sequence. As we mentioned in subsection 
IIIDl in the scenario at hand we can obtain the transi- 
tion from the matter-dominated era to the accelerated 
one. However, note that the dynamical system analy- 
sis can only give analytical results relating to the late- 
time states of the universe. The precise evolution of 
the universe towards such late-time attractors depends 
on the initial conditions, and it can only be investigated 
through a detailed numerical elaboration similar to the 
partial one that we performed in order to produce the 
aforementioned three figures. Thus, by suitably deter- 
mining the initial conditions, one can obtain universes 
that start with a de-Sitter expansion, then transit to the 
matter-dominated era, and finally falling into the late 
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FIG. 3: Invariant set x — 0, y — in the case of radiation 
(w — 1/3, n = 2). The shaded region corresponds to the 
unphysical portion of the phase plane. P^ {P2) is the local 



future (past) attractor, while B^ 
points. 



, B- and P{' , P^ are saddle 



time accelerating solution. The detailed examination of 
such behaviors, along with the investigation of the basins 
of attraction of the various evolutions in terms of the ini- 
tial conditions, and the estimation of the measures of the 
corresponding sub-spaces of the phase space, lie outside 
of the present work and are left for a future project. 

We close this section by discussing the cosmological 
evolution in the special case where Ei^ = Q an K = 
0, since in this case the Kantowsky-Sachs metric (in a 
comoving frame) 



ds" 



1) 



D' 



-dr' 



1 



D^{E^ 



rdr^ 



1 



3D'^K 



[di}^ +sm'^ if dip"^], (68) 



is singular. Although these points are unphysical their 
neighbouring solutions could have physical meaning, 
which can be extracted by obtaining first-order evolu- 
tion rates valid in an small neighborhood of the criti- 
cal points of the system. In particular, we first eval- 
uate the perturbation-matrix H at the critical point of 
interest, we diagonalize it and we obtain orders of mag- 
nitude for linear combinations of the vector components 
[Ei^jQ^'EjX^y) . Thus, the desired expansion can be 
obtained by taking the inverse linear transformation, 
and finally we preserve only the leading order terms 
as r — > — 00 or as r ^- 00. Although this procedure 
is straightforward in the case of dust matter (w = 0, 
n = 3/2) and radiation {w = 1/3, n = 2) we do not 



present it explicitly since we desire to remain in the gen- 
eral case of -Ei^ 7^ and K ^ 0. 



V. CONCLUSIONS 

In this work we constructed general anisotropic cosmo- 
logical scenarios where the gravitational sector belongs to 
the extended f{R) type, and we focused on Kantowski- 
Sachs geometries in the case of -R"-gravity. We performed 
a detailed phase-space analysis, extracting the late-time 
solutions, and we connected the mathematical results 
with physical behaviors and observables. 

As we saw, the universe at late times can result to a 
state of accelerating expansion, and additionally, for a 
particular n-range (2 < n < 3) it exhibits phantom be- 
havior. Additionally, the universe has been isotropized, 
independently of the anisotropy degree of the initial con- 
ditions, and it asymptotically becomes flat. The fact that 
such features are in agreement with observations [ij, |l6| 
is a significant advantage of the model. Moreover, in the 
case of radiation {n = 2, w = 1/3) the aforementioned 
stable solution corresponds to a de-Sitter expansion, and 
it can also describe the infiationary epoch of the universe. 

Note that at first sight the above behavior could be 
ascribed to the cosmic no- hair theorem [56[ , which states 
that a solution of the cosmological equations, with a pos- 
itive cosmological constant and under the perfect-fluid 
assumption for matter, converges to the de Sitter solu- 
tion at late times. However, we mention that such a 
theorem holds for matter-fiuids less stiff than radiation, 
but more importantly it has been elaborated for General 
Relativity [57|, without a robust extension to higher or- 
der gravitational theories |5al • In our work we extracted 
our results without relying at all on the cosmic no-hair 
theorem, which is a significant advantage of the analysis. 

Apart from the above behavior, in the scenario at hand 
the universe has a large probability to remain in a phase 
of (isotropic or anisotropic) decelerating expansion for a 
long time, before it will be attracted by the above global 
attractor at late times, and this acts as an additional 
advantage of the model, since it is in agreement with 
the observed cosmological behavior. However, the pre- 
cise duration of such a transient phase, unlike the at- 
tractor behavior, does depend on the initial conditions 
of the universe, which have to be suitable determined in 
order to lead to a deceleration-duration of the order of 
10^-10^*^ years, before the universe pass to the accelerat- 
ing phase. Such an analysis can only arise through an 
explicit numerical elaboration, that is beyond the ana- 
lytical, dynamical-system treatment which is where the 
present work focuses. 

The Kantowski-Sachs anisotropic i?"-gravity can also 
lead to contracting solutions, either accelerating or decel- 
erating, which are not globally stable. Thus, the universe 
can remain near these states for a long time, before the 
dynamics remove it towards the above expanding, accel- 
erating, late-time attractors. The duration of these tran- 
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sient phases depends on the specific initial conditions. 

One of the most interesting behaviors is the possibility 
of the realization of the transition between expanding and 
contracting solutions during the evolution. That is, the 
scenario at hand can exhibit the cosniological bounce or 
turnaround. Additionally, there can also appear an eter- 
nal transition between expanding and contracting phases, 
that is we can obtain cyclic cosmology. These features 
can be of great significance for cosmology, since they are 
desirable in order for a model to be free of past or future 
singularities. 

Before closing, let us make some comments concerning 
the use of observational data in order to constrain the 
present scenario. In particular, equation (|3T1l . along with 
the i?"-ansatz, relation (p4l) . and the definitions of the 
density parameters (l^5|) - (l^51) . can be straightforwardly 
written in the form used for observational fitting [59| . 
Thus, one can use observational data from Type la Su- 
pernovae (SNIa), Baryon Acoustic Oscillations (BAO), 
and Cosmic Microwave Background (CMB), along with 
requirements of Big Bang Nucleosynthesis (BBN) , to con- 
strain the model parameter n. Additionally, one can con- 
strain the initial allowed anisotropy value cr+ (or fio-), 
through its present value fio-o and the shear evolution 
([50]). Such a procedure is necessary for every cosmolog- 
ical paradigm, and can significantly enlighten the sce- 
nario at hand. However, since in the present work we 
focused on the dynamical features, the detailed observa- 
tional elaboration is left for a separate project. 

In summary, anisotropic _R"-gravity has a very rich cos- 
niological behavior, and a large variety of evolutions and 
late-time solutions, compatible with observations. The 
much more complicated structure of anisotropic geome- 
tries leads to radically different implications comparing 
to the simple isotropic scenarios. These features indicate 
that anisotropic universes governed by modified gravity 
can be a candidate for the description of nature, and de- 
serve further investigation. 
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Appendix A: Eigenvalues of the perturbation matrix 
H for the critical points 



of the perturbation matrix H calculated at each critical 
point. We also provide the eigenvalues of the perturba- 
tion matrix at the special points PI and P| located at 
the invariant circles Cg. 

For the two critical points Ac the associated eigenval- 
ues read 



x2l 



f 2(9w2 + 12w+l) 2(3w-|-2)(6u; + l) 
y^ i{iw + l) '^*^ 3(3w + l) 

(u; + l)(9u;2 + 12u; + l) 1 
~^ ZwTl J ' 

with [x2] denoting multiplicity 2. 

For the critical points Sg the associated eigenvalues are 



2(3u; + l) (3m; -2) (3m; -HI), , 2(m;+1) 
-e-r^ 7^, -e — [x2J,-e 



3(u;-1) ' 3(m;-1) 



w — 1 



For the critical points P^ the associated eigenvalues 
read 



0, e(3M; + 1), 3e(3M> -I- 1), --e{w - l)(3u; + l)\, 



while for P| they write 



0, 3e(3M; + l)[x2], --{w - l)(3w + 1) 



For the critical points P| the associated eigenvalues 
are 

{0,-e(3M;-t-l), 

9m;2 + V3V3u; + 1\/9m;3 + 21m;2 -f^ 31u; + 3 - 1 
~^ 2(3m;-1) ' 

9m;2 - ^/3^/3u; + 1V9m;3 + 21m;2 + ilw + 3 - 1 1 
~^ 2(3m;-1) J ' 

and thus there exists a l-dimensional center manifold 
tangent to the line 



X = -xQ 3m; + 1), y = e , ^-.w^ y-, 

2 iw -\- 1j(3m; — 1) 

S = -iQ(3M;-l),ge[-2,2] 



The system pSjl - lH^l) admits ten isolated critical points 
presented in Table |TT1 Here we provide the eigenvalues 



For the critical points PI the associated eigenvalues 
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6(3w + 1) (9^2 + 3u; + 1) 3(3w + 1) {9w^ + 3w + V9w^ + 3w + W45w^ + 51w + 5 + l) 



63u;2 + 30u; + 7 



63w2 + 30w + 7 



3(w + l)(3u; + 1) (9^2 + I2w + l) 3(3u; + 1) (9^2 + 3w - \/9w^ + 3w + lV45u;2 + 51u; + 5 + l) 
-e t:t: — ?; r:^ = -,—e ^^ 



63w2 + 30w + 7 



63u;2 + 30w + 7 



Finally, for the critical points P^ the associated eigen- 
values read 



0, 0, -e '-- '-, 6e{w + 1 

3z« — 1 



Appendix B: The direction Ei^ 



In subsection IIIICI we have investigated the system 
(|i5)l - (|i5)l . leaving outside the decoupled equation (150)1 . 
However, this equation provides information of how 
Kantowsky-Sachs f{R) models are related to more gen- 
eral anisotropic geometries. For instance from ()50p we 
see that Ei^ = is an invariant set of (H51) - ([5D|) . This set 
is closely related with the so-called "silent boundary" one 
p3(]| - l62| . ^ Note that, as current investigations suggest, 
in Go-cosmologies, that is in cosmologies where there is 
no symmetry, both past and future attractors belong to 
the silent boundary [6l|, |6J] . 

In the present work we do not focus in how Kantowsky- 
Sachs f{R) are related to more general anisotropic mod- 
els, and thus, we do not elaborate completely equation 
((50)) . However, we can still obtain the corresponding sig- 
nificant physical information, since the stability along the 
Ei^ direction is determined by calculating the sign of 
d {Ei^) /dEi^ , which coincides with the eigenvalue asso- 
ciated to the eigenvector Ei^. Thus, if it is negative then 
the small perturbations in the Ei ^ direction decay, while 
if it is positive they enhance. 

From the sign oi d (E,^)' /dE.^U = - ^'ttlTij'^^ 
it follows that A+ {A- ) is stable (unstable) to small per- 
turbations in the Ei^ direction provided w+ < w < 1, 
otherwise it is unstable (stable). 

^ The concept of "silent boundary" was introduced for inhomo- 
gcncous cosmologies and it is defined by the condition e\/j3 = 
where P = H + a^ , with the /3-gradient being equal to zero. It 
corresponds to the divergences of H + cr_|_. In particular, as the 
initial singularity is approached, the radial part of the metric 
tends to zero and the null cone becomes a timeline (see section 



From the sign of d (Ei^)' /dEi^B^ = -f^^J^ ^^ fol- 
lows that S+ (S-) is unstable (stable) to small pertur- 
bations in the Ei^ direction. 

From the signs of d (Ei^)' /dEi^\p^ = 2e(3w-f 1) it 
follows that Pj^ (^1^) is stable (unstable) to small per- 
turbations in the Ei^ direction provided — | < w < 1, 

otherwise it is unstable (stable). Since d (^Ei^^ /dEi^ 
vanishes at P2 and at P|, it follows that they are neu- 
trally stable to small perturbation in the Ei^ direction. 

By analyzing the sign of d (Ei^) /dEi^\pc = 

"'eaJii+aou^+T' — "^^ deduce that P4+ (P;^) is sta- 
ble (unstable) to perturbations in the Ei^ direction if 
w+ < w < 1. 

From the sign of d {Ei^)' /dEi^\p. 
(3.+r)((3.-i).+.V3^) ^ .^ ^^^^^^^ ^^^^ ^^; ^^_^^ 

whenever exists, is always unstable (stable) to small 
perturbations in the Ei ^ direction. 
Finally, note that since 



d{E,^y 

dEi^ 



\c, =2(n-l)[e-|-cos(w)]. 



then for 1 < n < 3, C+ is always unstable (except for 
M = tt) to small perturbations in the Ei^ direction, while 
C_ is always stable (except for u G {0, 2tt}). 

Lastly, it is interesting to mention that all the critical 
points of the reduced system, except PI 3, belong to the 
"silent boundary" in the full phase space, that is each one 
has a representative in ^ x R with Ei^ — 0. In particular, 
the critical points Ae, Be correspond to isotropic silent 
singularities of the full five-dimensional phase space and 
P^ correspond to an anisotropic one. 

5.3 of 1631 for a clear physical discussion). In our analysis we do 
not use /3-normalization, however, as E^ = is approached, a 
similar physical behavior is attained (that is the radial part of 
the metric tends to zero and the local null cone collapses onto 
the timeline). 
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